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We consider the entanglement properties of free fermions in one dimension and review an 
approach which relates the problem to the solution of a certain differential equation. The 
single-particle eigenfunctions of the entanglement Hamiltonian are then seen to be spheroidal 
functions or generalizations of them. The analytical results for the eigenvalue spectrum agree 
with those obtained by other methods. In the continuum case, there are close connections 
to random matrix theory. 

I. INTRODUCTION 

The entanglement properties of many-particle quantum states have been the topic of many 

n 

studies in recent years [1]. This holds in particular for the ground state of free fermionic systems. 
These are Slater determinants, and in this case the reduced density matrix (RDM) for some portion 
of the total system has the form 

P=\^ (1) 

where % is again a free-particle Hamiltonian, see Q|. Its single-particle eigenfunctions (p^ and the 
corresponding eigenvalues et can be determined either from the one- par ticle correlation functions 
[3-^] or from the overlap of the occupied states in the subsystem \\. lid- 12]. Their properties for 
the most studied case, namely a segment in a chain of free fermions, are relatively well-known. 
Thus the low-lying \Ek\ vary roughly linearly with k and the slope is proportional to 1/lnL if the 
subsystem has length L and InL is large. This leads to a logarithmic variation S = 1/3 In L of 
the entanglement entropy S which is characteristic for critical systems and follows from conformal 
invariance, see [3]. The law was first obtained from an asymptotic analysis of the correlation 



matrix using the Fisher-Hartwig conjecture 
in the same way 



15, 



14l ]. and further subleading terms have been derived 
1.0], I ! ' (ugcnl'uuot ions (fk are largest near the boundaries for small \ek\ and 
concentrated in the interior of the subsystem for large \£k\. This can be seen easily from numerical 
calculations. For a half-filled system, there are also analytical expressions in the continuum limit 



1. 



9]. Altogether, a picture emerges, where the (entanglement) Hamiltonian H has both bulk and 
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boundary-like states with the latter dominating the entanglement. 

Given the relatively simple form of the correlations (or the overlap), one wonders whether a 
complete analytical treatment of the problem is possible. This is, in fact, the case and it turns 
out that it was done a long time ago by Slepian in the analysis of time and band limited signals 



17, ll8|]. A good account of the background can be found in his 1982 John von Neumann lecture 
191 ] . The work has received a large number of citations in very different areas. The analogy to the 
entanglement problem, where one deals with limited regions in momentum space (the Fermi sea) 
and in real space (the subsystem), was pointed out already by Gioev and Klich [8]. However, while 
they derived a formula for S in higher dimensions, we focus here on the eigenvalue problem. 

It turns out that the eigenfunctions (fk, or their Fourier transforms, are the (prolate) spheroidal 
functions which appear if one separates the Helmholtz equation in elliptical coordinates, or general- 
izations of them. They enter in the continuum case because the integral kernel of the entanglement 
problem commutes with the corresponding differential operator. Those concentrated in the interior 
of the subsystem are quite familiar objects, namely oscillator functions. They become more com- 
plicated once they touch the boundary. In any case, the machinery of solving differential equations 
can be used to obtain their form and also the eigenvalues e^. In particular, an asymptotic analysis 
is possible and gives the low-lying ones for the case of a large subsystem. This leads to exactly 
the same density of states as found in the Fisher-Hartwig approach Q, l^ . If either momenta 
or positions become discrete, the situation is still similar. Working with the other, continuous 
quantity, one has a commuting differential operator and can discuss its eigenfunctions. One should 
mention that a lot of this material also appears in random matrix theory, because in the Gaussian 



unitary ensemb 
interest, see 



e the same integral kernel enters and the eigenvalues Ek determine the functions of 



211 ] . The connection to this field was pointed out by Keating and Mezzadri 22], [23 ] 
and also used in [2], but again with the focus on the entropy. 

The purpose of this paper is to draw attention to the results described above, because they 
complement the usual approaches and put the entanglement problem into a broader context. Our 
own contribution is mainly the collection and the presentation, including a number of figures. In 
section 2 we give a brief general outline of the determination of <pk via correlation and overlap 
matrices. In section 3 we discuss the case of an infinite continuum system, where the simple 
sine kernel and the usual spheroidal functions appear, of which we show some examples. Section 
4 treats the case of a finite continuum system, where the correlations still lead to an integral 
equation. Section 5 deals with an infinite lattice where the same equation appears in momentum 
space. In this case we also discuss the matrix which commutes with the correlation matrix, as well 
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as the dispersion relation of the sj.. Finally, section 6 contains remarks on semi-infinite systems 
and higher dimensions and a brief conclusion is given in Section 7. 

II. CORRELATION MATRICES AND OVERLAP MATRICES 

The structure of the problem is best seen from the formulae for a general lattice system with a 
discrete set of orthonormal single-particle functions <& g (n). A many-particle state \F) in which a 
set F of these is occupied (the Fermi sea) then is 

\F}=Uc\\0) (2) 

q&F 

where c\ are creation operators and |0) is the vacuum. With 

q 

the correlation matrix in state \F) becomes 

C mn = (F\clc n \F) = *;M* ff (n) (4) 

q&F 

and the RDM in a subsystem S is determined by the eigenvalue problem 

^2 Cij ipk{j) =Ck<Pk(i), i G S (5) 
jes 

If S consists of L sites, ([5]) gives L single-particle eigenfunctions ip^. The eigenvalues Cfc with 
< (k < 1 are their non-integer occupation numbers and related to the 8k via 

£h = In — — — or Ck = (6) 

These quantities can also be obtained by constructing the Schmidt decomposition of \F) directly Q]. 
In this case, one forms new orthonormal functions from the occupied Q q which are orthogonal 
also in the subsystem (and in the remainder R). This is done by calculating the overlap matrix 

and solving the eigenvalue problem 

^ A qq > ip k (q') = (k<fk(q), q^F (8) 

q'&F 
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Then 

tf fc (i) = X^Gz)**® (9) 

satisfies when restricted to S. Similarly, its restriction to R satisfies the analogue of ([5]) with 
eigenvalue (1 — Moreover, it has norm in S and norm 1 — ^ in R. Therefore, if ipk is 
normalized in S, one has 

<p k ({) = -^^k(i), i G 5 (10) 

which gives the pair of relations 



^ q eF v Cfc ig5 

The normalization properties allow to obtain the eigenvalues solely from the Vljfc if these can be 
found by some other means (as will be the case below). If one arranges the Cfc in decreasing order, 
the corresponding ^ are less and less concentrated in S. For L sites in the subsystem, ([8]) can 
have only L non-trivial solutions. If the particle number, i.e. the number of states in F, exceeds 
L, the remaining ^k have no components in the subsystem and = 0. 

In the following we apply these formulae to the case of plane waves, 3> 9 (n) ~ e ign or & q (x) ~ e iqx . 
The two eigenvalue problems (|5|) and ([8]) then correspond to working in real space and in momentum 
space, respectively. Moreover, the state \F) will be a normal Fermi sea (\q\ < qp) and the subsystem 
S a single segment. The two functions (fk{i) and <fk(q) are then related by a Fourier transform 
limited to a window in q space and in real space. 

III. INFINITE CONTINUOUS SYSTEM 

We now consider free fermions on an infinite line with all momentum states \q\ < qF filled. The 
system then has average density n = qF/it and the correlation matrix is 

C{x - x>) = f 9F * = ™QF(x-f) (12) 

J-q F 2vr vr(x - x') 

Choosing the subsystem 5 as the segment — 1/2 < x < £/2, the eigenvalue problem © becomes 

i-e/2 

dx'C(x - x')cpk(x') = Ckfk(x) (13) 

1/2 

or in reduced variables y = 2x/t with (fk(ly/2) = ipk(y) 

1 dy'K(y-y')My')=CkMy) (14) 
l 
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with the sine kernel 



, smcjy - y 

K{y -y) = -, jr- , c = q F £ 2 (15) 

vr(y - y') 



Similarly, the overlap matrix is 



A(q - q') = I"'" p. = ^ = 4) H 2 (16) 

' /2 2vr 7r(g- g ') V ' 



and with p = q/qF the eigenvalue problem ([8]) becomes 

l 

dp' K(p - p')ip k {p) = Ck ^k{p) (17) 

-l 

Thus the equations in real space and in momentum space are identical. In the following, we will 
work in real space. 

The sine kernel is the "square" of another, even simpler symmetric kernel K(y, y') 



K(y-y') = j\zK(y,z)K*(z,y'), K{y,z) = (18) 



Therefore, the eigenvalues Ck of the sine kernel follow from the eigenvalues \Xk of K via Ck = |^fc| 2 - 
Moreover, K commutes with the second-order differential operator 

D -Ty {1 ~K + ^ <19) 
provided the functions one operates on are regular at y = dbl. Therefore the (real) eigenfunctions 
of K and K can be obtained from 

£>V = H (20) 

This equation has continuous and bounded solutions for all y only for a discrete set of values 9 = 9k- 
The correspondin g ij^ u are the prolate spheroidal wave functions of the first kind and order m = 0. 
In the notation of 25f] these are called angular functions S$k f° r y 2 < 1 an d radial functions i?ofc f° r 
y 2 > 1 (which fits with our notations of subsystem and remainder) . Both are normalized separately 
and connected by joining formulae. Their properties have been investigated in great detail, see e.g. 



26H28t] . Studying the solution of (|20p for all y gives the quantity ^ of (|10p . Of particular interest 



is the case where the subsystem contains many particles. Since the average number is 



N = in = -c (21) 

7T 



this corresponds to large c. 
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A qualitative picture of the spheroidal functions can be obtained by observing that for small y 
the operator D reduces to 



dy 2 



D = -— I +cY (22) 



and thus to the Hamiltonian of the harmonic oscillator with frequency u = c. The oscillator 
functions 

u k (y) =A k e-^y 2 H k (^y) (23) 

with the Hermite polynomials H k have a spatial extension y k ~ sj (2k + l)/c (using the classical 
turning point) and therefore "fit" into the subsystem as long as 2k + 1 < c. Under this condition 
the exact spheroidal functions resemble them closely, in particular those concentrated fully in the 
interior. A picture gallery is shown in Fig. [IJ where we plotted the first 15 spheroidal and oscillator 
functions for c = 5tt (N = 10). Both are normalized to 2/ (2k + 1) in (—1, 1). The function Sok has 
exactly k zeros there and is alternatingly even and odd. One sees that up to k ~ 7 the differences 
are rather small. As k increases further, the oscillator functions deviate more and have zeros 
outside the subsystem, while the oscillation of So*; becomes faster near the boundaries. For large 
c, it is logarithmic up to a distance from ±1, i.e. besides y the variable ln[(l + y)/(l — y)] 

enters, as also found in a continuum approximation for a half-filled lattice system [9]. The nature 
of the eigenfunctions is reflected directly in the occupation numbers Qk which are close to 1 for 
small k and close to zero for large k. As discussed below, the transition takes place at k ~ N. 

If one continues tp^ to y 2 > 1, one finds the eigenfunctions in the remainder R corresponding 
to the eigenvalue 1 — £fc- These are a kind of mirror image of those in S. For small k, they are 
oscillating rapidly near the boundary but have very small amplitudes (the phenomenon has been 



termed superoscillation 291]). while for large k the behaviour is smooth there and only later an 
oscillation sets in. This is illustrated in Fig. [2] where the i?ofc are shown for a number of k values 
up to k = 40. Note that due to the normalization conventions the functions S'ofc and i?ofc do not 
match at y = 1 The form for \y\ S> 1 can be obtained from D by writing ip(y) = x(v)lv which 
leads to the operator 

* _2 



and gives approximately wave-like solutions x(y) ~ cos(cy + a) = cos(qpx + a). This can be seen 
clearly in the figure. 
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FIG. 1: The first 15 spheroidal functions Sofc(c, y) for c = 5tt (full lines, red) and, for comparison, the 
oscillator functions Uk(y) (dotted lines, blue). 

Turning to the eigenvalues (k (usually called in the literature), they can be obtained from 
the relation Q 

J dye^'Sokfay) = 2i n R 0k (c,l)S 0k (c, z) (25) 
which is, up to a factor, the eigenvalue equation of the integral operator K. This gives the general 
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formula 

Ck = W\ 2 = -[Rok(c,l)} 2 (26) 

7T 

which can be evaluated numerically. To obtain analytical expressions for large c, Slepian considered 
the solutions of (|20p in the various regions of y and connected them in the domains of overlap [l^] • 
While he considered small and large k separately, des Cloizeaux and Mehta using WKB formulae 
could later obtain a single expression 



30j] and Landau and Widom gave a mathematical proof j^j. 



The procedure is rather tedious, but in the end a relatively simple result for emerges if k ~ N, 
c.f. (1.37), (1.38) in [l?]]. Namely, it is given by the root of smallest absolute value of the equation 

c+| H^c)-^) = \{k-\) (27) 

where <p(z) = argr(l/2 + iz) and T(z) is the gamma function. This has solutions <C 1 which 
are obtained by expanding (p(z) 

g(ln(4c)-^(0)) = |(fc-l/2-iV) (28) 

or, with c inserted, 

**(k-l/2-N) 
k ]n(2wN) - V>(l/2) 1 1 

where ip(l/2) = —7 — 2 In 2 = —1.963..., 7 is Euler's constant and ip(z) denotes the digamma 
function. Thus et goes through zero through 1/2) for k ~ N 3> 1 and the slope vanishes 
logarithmically with c or N. A similar formula was found in [l| for a half-filled lattice system. For 
larger one really has to solve ()27|) . It is simpler, however, to view it as a relation k = k{e) which 
leads to a density of states n{e) = dk/de given by 

' In I \r) - :'l 



n e = — 

7T Z 



l n(4c) - y/( A)l (30) 



where explicitly 



^) = ^[^ + «)+^-w)] (31) 

This is the result found and plotted in [20] and implicit already in [l^ |. Inserting it into the 
expression for the entanglement entropy, one finds the leading logarithmic term and the first 
correction. We will come back to it in Section 5. 
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FIG. 2: Spheroidal functions Rok(c, y) for c = 57r, 1 < y < 3 and various k. 



IV. FINITE CONTINUOUS SYSTEM 



Consider now a finite system in the form of a ring with circumference L. The momenta are 
quantized according to q = 2nn/L with integer n and the Fermi momentum is written as qF = 
2nm/L. This gives a total particle number N = 2m + 1 and a correlation matrix 



C(x - x') 



1 



E 



-i2-K(x-x')n/L _ 1 sin (Nn/L)(x - X 1 ) 



L ^— ' L sin (ir/L)(x — x') 

n=—m \ i / \ / 



(32) 



With the same subsystem as before, — £/2 < x < £/2, the eigenvalue equation (|13p retains its form. 
However, here it is preferable to define the reduced variable as z = x/L such that \z\ < 1/2. Then 
the equation becomes 



w 



dz'K'(z-z')Mz') = CkMz) 



(33) 



with W = £/2L < 1/2 and the modified sine kernel 

simrN(z — z') 



K\z - z') 



sin7r(z — z') 



(34) 



10 



The overlap matrix is 



{q q) ~J-mL e ~ (q-q')L/2 



-1/2 

but since the momenta are discrete, one can write it as 



(35) 



sin (nl/L)(n-ri) 

A nn ' = J (36) 

■K{n — n) 

and the eigenvalue equation is a genuine matrix equation with a Toeplitz matrix. 

The problem is now in exactly the form as studied by Slepian in 1978 [la ]. In particular, (|33j) 
is his equation (10) and the ipk are bis functions U/-. There are only N of them with non-zero 
eigenvalue (corresponding to the N eigenvalues of the overlap matrix) and he calls them discrete 
prolate spheroidal wave functions. The treatment of the infinite case can be repeated, because 
there is again a commuting differential operator, namely 

D' = --^^-(cos2ttz -cos2ttWO-^ - \{N 2 - l)cos2vrz (37) 
4ir 2 dz dz 4 



One sees that, in comparison with (j 1 9 [> . the powers have been replaced by cosine functions, which 
reflect the ring geometry. The first term is a kinetic energy with different sign in S and R, as in 
the infinite case. The second term is a potential with a minimum at the origin. In the limit of 
small z and W one can expand the cosine functions and recovers the operator D of the infinite 
case. The finite size is then no longer visible. The same holds for the integral operator K' which 
reduces to K. The TV" lowest eigenfunctions of the differential operator give the solution both in 
the subsystem and in the remainder, and determining the respective norms allows to calculate the 
eigenvalues £fc- Pictures of the functions for N = 4, 5 and W = 0.2,0.4 are shown in [la] and one 
sees the similarity to the oscillator functions in particular for the lowest state k = 0. 

The analysis of the eigenfunctions of D' for large N leads to formulae for which are very 



similar to those of the infinite system, see equs. (53), (60) and (61) in 18]. Instead of (f27|) one 
finds the equation 

itNW + ^ ln(2iV sin 2irW) - <p(^) = ^-(k-\) (38) 

Z7T Z7T 2 2 

and in the linear region this gives after inserting W 

nHk-l/2-N) 
k \n{2Nsm{ni/L))-i){l/2) K ' 

For L — > oo, keeping N/L = n constant, this reduces to (|29p . 

The existence of the commuting operator D' also has consequences for the other eigenvalue 

problem. This will be discussed in the following section. 
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V. INFINITE LATTICE SYSTEM 

In the entanglement investigations, this is the most studied case. Let the system be the limit 
of a ring with a total of M sites. The correlation matrix is 

C nn ,= T ^e-^ = Sm f n -^ (40) 
J-q F 2vr ir(n - ri) 

which corresponds to (|12p but with discrete sites. We will choose the subsystem here as the L sites 
n = 1, 2, ...L. Then the eigenvalue equation is 

L 

^ Cnn' Pk(jl') = ( k <Pk(n) (41) 
n'=l 

and for the overlap matrix one finds 

L 



n=l 

?(L+l)/2 ( 1 Sm(g-g0V2 ^ p ig'(L + l)/ 2 

VM shx{q-q')/2 



e iq'(L+l)/2 ^ 



The exponential factors reflect the location of the subsystem, but do not affect the spectrum. The 
eigenvalue equation can be written 

r9F d -f sin( r g/ iw ? 2 *>w = & ^ = e " (L+i)/2 ^ w 

-q F 2?r sm[q - q')/2 

With p = q/2-K, W = qF/2it = n/2 and <p k (2Trp) = tp(p) this leads to 

-w 

dp' K'(p - p')ip k (p') = (k A (p) (44) 

w 

re kernel K'. 



which is ()33p with the substitution N — > L in t 
Therefore, as noted in previous work 10J, 



111 ], the infinite lattice and the finite continuum 



problem are actually the same with the role of positions and momenta, i.e. of correlation and 
overlap matrices, interchanged. Thus the discrete spheroidal functions appear here in momentum 
space, and the same holds for the commuting differential operator. As mentioned above, this 
operator has a consequence for the other eigenvalue problem ([41 p . From ([lip 

<p k {q) = -±_ e ^L + i)/2 1 £ Mn)e ^n (45) 
VCk VM ^ 

and inserting this into the eigenvalue equation D'ip = 0'(p one finds a three-term recursion relation 
for the (pk{n) which involves only the sites n— 1, n and n + 1. In other words, the (fk( n ) are also 
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the eigenfunctions of a tridiagonal matrix T which therefore commutes with C (in the subsystem) . 
Thus the relation [K',D'] = has an equivalent [C, T] = in real space. 
Writing T in the form 

/ di h \ 

ti d 2 t 2 

I t 2 d 3 t 3 (46) 



tL-i d L J 



the matrix elements are 



t n = -^n(L - n), d n = ( 



L + l 



n) 2 cos qF 



(47) 



Both coefficients vary parabolically around the center of the system. Viewed as a hopping Hamil- 
tonian, T describes an inhomogeneous system where the hopping is largest in the middle. The 
potential energy is also largest there if qF > vr/2, but smallest if qp < vr/2. For a half-filled system, 
qp = vr/2, it vanishes. 

This matrix was later rediscovered in [9] on the basis that the Hamiltonian in ([1]) when calculated 
numerically shows a very similar structure. Note that T is only determined up to a multiplicative 
factor and an additive constant. The matrix elements given in [9J differ from (|47|) by the factor 2/L 2 
and a constant in d n . Using T, one can obtain the eigenfunctions (fik numerically for much larger 
systems than via the correlation matrix. They resemble the spheroidal functions of the continuum 
and were termed discrete prolate spheroidal sequences by Slepian. In Fig. [3] we show some of them 
for a system of L = 50 sites and qp = vr/2. For qp = 0, vr they are given by Hahn polynomials and 
the eigenvalues of T are very simple [32]. In the limit qp — > and L — > oo keeping qpL constant, 
the problem becomes continuous in y = n/L and one comes back to the usual spheroidal functions. 
Expanding the quantities in T<p = 9'ip, one can also derive the differential operator D directly. 

The eigenvalues e& follow directly from the formulae in section 4 with N — >• L and 2nW = qF. 
Then ([39]) becomes 

TT 2 (k-l/2- N) 



(48) 



ln(2Lsin q F )-i/>(l/2) 

For a half-filled system (qp = tt/2) this gives the asymptotic result found in Q], whereas for qp — > 
and L — > oo it reduces to (|29p . The density of states n(e) which decreases with e is clearly seen in 
numerical calculations of the 2\. This is illustrated in Fig. [U where the eigenvalues obtained 
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FIG. 3: The first 15 eigenfunctions ifk for a subsystem of L = 50 sites in a half-filled infinite chain. 

directly from (|4ip are shown together with the theoretical formula for k = k(e), viewing A: as a 
continuous variable. The curves are not exactly linear, but show an upward (downward) bend on 
the positive (negative) side. One sees that for L = 20 there are small differences between both 
results, but for the larger L there is an excellent agreement. One also sees that even there the 
linear region is still quite small. 
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-15 -10 -5 5 10 15 

k- (L+l)/2 

FIG. 4: Dispersion relation of the for a half-filled infinite lattice and three different sizes of the subsystem. 
Symbols: numerical data, lines: analytical result for k(s). 

VI. FURTHER ASPECTS 



So far, we have considered only systems without boundaries. However, the results cover also 
the case of a semi-infinite geometry when the subsystem is located next to the boundary. The & q 
are then sine functions and, in the continuum case, lead to the correlation matrix 

C a [x, at) = C(x-x') -C(x + x'), x,x'>0 (49) 

where C(x — x') is the result (|12p. But if the subsystem is the interval < x < £/2, the eigenvalue 
problem 

/ dx' C s (x - x')ifk(x') = Ck ¥k{x) (50) 
J o 

is solved by the odd eigenfunctions of (|13|) . Thus one can take over the previous results, confining 
oneself to the odd fe-values. This means in particular that the spacing of the low-lying doubles 
and the density of states halves. As a result, the entanglement entropy is also halved. 

The argument also holds for the lattice case, where C s has the form f)49|) with x, x' — > n, n' 
and n,n' > 1. However, the solutions of the analogue of (j50[) then vanish at the point n = 0. 
Therefore, one has to compare with an infinite lattice problem, where this point is the center of the 
subsystem. This means, that it has to consist of 2L + 1 sites (an odd number), if in the semi-infinite 
case the subsystem has L sites. The connection between infinite and semi-infinite system is also 
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seen at the level of the commuting matrices. There exists a tridiagonal matrix T s with the property 
[C S ,T S ] = and it is just the right half of the matrix T 

Returning to closed systems, one should mention that a finite ring of M lattice sites filled with 
N = 2m + 1 particles leads to the correlation matrix (|32p with x — > n and L — > M. For a segment 



331 ] which is a generalization of T 



of L sites, one can then also find a commuting tridiagonal matrix 
and leads to another discrete version of the spheroidal functions. The parabolic form of the t n and 
d n is then replaced by trigonometric expressions. For example t n ~ sin(7rn/M) sin(7r(L — n)/M), 
and for M — > oo one reobtains (|47p . 

Another generalization concerns higher dimensions. The case of an infinite continuum with 
spherical Fermi surface and a spherical subsystem has also been studied in detail [34]. Due to 
the rotational symmetry, the eigenfunctions are products of an angular and a radial factor. For 
the latter, one finds an integral equation with a Bessel function as kernel. Again a commuting 
differential operator exists and is closely related to D. For d = 2 and with y — > r 

D r = D+ m2 ~ 1/4 , r > (51) 

where m = 0,±1,±2... is the azimuthal quantum number. Physically, the additional term can be 
viewed as a centrifugal potential. As a result, the eigenfunctions have to vanish at the origin. They 
were called generalized spheroidal functions and some are shown in [3^]. For the eigenvalues e^ m 
for fixed m one finds very similar results as in one dimension. 



VII. CONCLUSION 



We have considered the ground state of free fermions and discussed an approach which deter- 
mines the reduced density matrix directly via its single-particle eigenfunctions. It circumvents the 
original matrices or integral kernels and works instead with differential operators. In the continuum 
case, this leads to spheroidal functions which are well known in mathematics and appear in various 
other areas. Physically, they have the character of either bulk or boundary-like states, and it is 
amusing that, in the first instance, they are close to harmonic oscillator functions. In hindsight, 
however, one could have guessed that from the appearance of their lattice counterparts. In the 
lattice case, the asymptotic results for the single-particle eigenvalues bridge a gap between the 
numerical calculations and the Fisher-Hartwig determinantal formulae which only give a density 
of states. We have demonstrated in Fig. [J] how well the results match. 

On the technical level, the approach is rather tedious and inferior to Fisher-Hartwig calculations. 
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However, the concept of a commuting operator is interesting and reminiscent of transfer matrix 
problems in integrable classical lattice models. It is also interesting that this feature carries over to 
higher dimensions for spherical subsystems. As to the one-dimensional problem, one has now an 
almost complete overview over the properties of the entanglement Hamiltonian H, but an explicit 
form is still lacking. 
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